Performance of Wang-Landau algorithm in lattice model of liquid crystals 
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We present a study on the performance of Wang-Landau algorithm in a lattice model of liquid 
crystals which is a continuous lattice spin model. We propose a novel method of the spin update 
scheme in a continuous lattice spin model. The proposed scheme reduces the autocorrelation time 
of the simulation and results in faster convergence. 
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I. INTRODUCTION 

The Wang-Landau (WL) algorithm [H, introduced in 
2001, has received much attention and has been applied 
to a wide range of problems • In most of these inves- 

tigations, the authors have applied the WL algorithm to 
systems with discrete energy levels. However, relatively 
fewer papers have so far appeared on lattice models with 
continuous energy spectrum 0, fT3l - [l7| . Techniques, in 
general, to improve the algorithm for different problems 
have also been proposed dS'^S^. The review [s^l illus- 
trates the versatile applications of the WL algorithm in 
protein folding, fluid simulations, systems with first order 
phase transitions and other systems with rough energy 
terrain. Some authors find its applications in performing 
numerical integration [13, 113 • 

The WL algorithm allows us to calculate the density of 
states (DOS) as a function of energy or the joint density 
of states (JDOS) as a function of energy and a second 
variable [ij]. For a macroscopic system, the DOS fl{Ei) 
(where i = 1,2, ■ ■ ■ ,n, n being the bin index) is a large 
number and it is convenient to work with its logarithm 
g{Ei) = In Q{Ei). Since the DOS is independent of tem- 
perature and contains complete information about the 
system, the task is to determine it as accurately as pos- 
sible. The next step involves the determination of parti- 
tion function Z{T) = e"''^ (/3 = 1/T, Boltz- 
mann constant has been set to unity) at any temperature 
(T) by the standard Boltzmann reweighting procedure. 
Once the partition function is known, the model is es- 
sentially "solved" since most thermodynamic quantities 
at any temperature can be calculated from it. The algo- 
rithm is implemented by performing an one-dimensional 
random walk that produces a "flat" histogram in the en- 
ergy space. For a continuous model, one needs to use a 
discretization scheme to divide the energy range of inter- 
est into a number of bins which label the macrostates of 
the system. In the WL algorithm, these macrostates are 
sampled with a probability which is proportional to the 
reciprocal of the current DOS. The estimate for the DOS 
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is improved at each step of the random walk using a care- 
fully controlled modification factor / to produce a result 
that converges to the true DOS quickly. A histogram 
record H{Ei) of all states visited is maintained through- 
out the simulation. When g{Ei) corresponding to a cer- 
tain macrostate is modified as g{Ei) — ^ g{Ei) + In/, the 
corresponding H{Ei) is modified as H{Ei) H[Ei) -\- 1. 
In the original proposal of WL algorithm, an iteration is 
said to be complete when the histogram satisfies a certain 
"flatness" condition. This means that H{Ei), for all val- 
ues of i, has attained 90% (or some other preset value) 
of the average histogram. In the following iteration, / 
is reduced in some fashion, the H{Ei)'s are reset to zero 
and the process is continued till In/ is as small as 10^^ or 
10~^. Since the history of the entire sampling process de- 
termines the DOS, the WL algorithm is non-Markovian 
besides being multicanonical in nature. 

In course of the random walk in a WL simulation, the 
fluctuations of energy histogram, for a given modification 
factor /, initially grows with time and then saturates to 
a certain value. Zhou and Bhatt [131 carried out a math- 
ematical analysis of the WL algorithm. They provided a 
proof of the convergence of the iterative procedure and 
have shown that the fluctuations in histogram, propor- 
tional to l/v^hT/ for a given /, cause statistical errors 
which can be reduced by averaging over multiple simula- 
tions. They have also shown that the correlation between 
adjacent records in the histogram introduces a system- 
atic error which is reduced at smaller /. The prediction 
m Ref. [271 has been numerically verified by different au- 
thors independently p/7,, ,29] . Although to obtain a flat 
histogram is the initial motivation behind the WL algo- 
rithm, Ref. [13 concluded that flatness is not a neces- 
sary criterion to achieve convergence and suggested that 
one should instead focus on the fluctuations of the his- 
togram rather than the "flatness" . They had shown that 
l/y/TnJ visits on each macroscopic state is enough to 
guarantee the convergence. In fact, fluctuations in the 
histogram is intrinsic to WL algorithm. These fluctua- 
tions lead to a statistical error in the DOS which scales 
as \/hi f, for a given /. The iterative WL algorithm par- 
tially reduces this statistical fluctuations by decreasing / 
monotonically. However Ref. [Hj clearly illustrates that 
^^^hSSB^l / is reduced to a very small value according to 
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the original prescription, the statistical error stops to de- 
crease at a certain point. In practice there always exists 
a systematic error in the simulation which is a function 
of / and the correlation between adjacent records in the 
histogram. Ref. ^] observed that this systematic error 
decreases when either / or the correlation decreases. In 
this context, we refer to the work of Morozov and Lin 
3l| who presented a study on the estimations of accu- 
racy and convergence of the Wang-Landau algorithm on 
a two level system with a significant efficiency improve- 
ment in [sHi. The WL algorithm compares r2(i?i) and 
^l{Ef), i.e, DOS before and after an attempted move, 
but it does not require Ei to be close to Ef. This is why 
Rcf. [13] suggested the use of cluster algorithms that 
allow "nonlocal" moves in the parameter space. The 
Ref. I'ightly pointed out that the update schemes 

for the underlying model certainly have an effect on the 
outcome. In the present paper we suggest a method for 
the spin update scheme of a lattice model with continu- 
ous energy spectrum, which reduces the autocorrelation 
time by an appreciable amount compared to the conven- 
tional spin update scheme. The suggested spin update 
method to obtain a less correlated configuration has also 
the advantage that this method is free from tuning any 
adjustable parameter. The method is described in Sec- 
tion [Till We also investigate the growth of the histogram 
fluctuations in the one-dimensional Lebwohl-Lasher (LL) 
model, described in Section [ill to check if the nature of 
the dependence of the maximum of the histogram fluctu- 
ations on the modification factor / is model independent 
or not. We mention in passing that Ref. 0] suggested 
the model-independent nature of the maximum of the 
histogram fiuctuations by performing simulations on two 
discrete Ising models and concluded that many more sim- 
ulations on different models are needed to confirm this 
universality nature. Ref. (l7| confirmed this universality 
behavior for two continuous lattice spin models with spin 
dimensionality two. We have found that for the present 
model (spin dimensionality three) , the fiuctuations in the 
energy histogram, after an initial increase, saturates to a 
value which is inversely proportional to v^hT/ and con- 
firm that this feature is generic to the WL algorithm. In 
the second part of the work, we have carried out the WL 
simulation with the proposed spin update scheme to es- 
timate the canonical averages of various thermodynamic 
quantities for lattices of reasonably large size where min- 
imum number of visits to each macrostate are l/\/\nJ. 
Results obtained from our simulation are compared with 
the exact results available for the model. 



The rest of the paper is arranged as follows. In Sec- 
tion ini we have described the model. The computa- 
tional techniques are discussed in Section IIIIl Section IIVI 
presents our results and discussions. Section |V] draws the 
conclusions. 



II. MODEL 

For the purpose of investigation, we have chosen an 
one-dimensional array of three-dimensional spins {d = 
1,^ = 3, where d is the space dimensionality and I is the 
spin dimensionality) interacting with nearest neighbors 
(nn) via a potential 



V^J - -P2(C0S e,j) 



(1) 



where P2 is the second Legendre polynomial and 9ij is 
the angle between the nearest neighbor spins i and j 
(the coupling constant in the interaction has been set 
to unity). The spins are three-dimensional and head- 
less, i.e, the system has the 0(3) as well as the local Z2 
symmetry, characteristic of a nematic liquid crystal. The 
model, known as the Lebwohl-Lasher (LL) model f36l. is 
the lattice version of the Maier-Saupe (MS) model |37| 
which describes a nematic liquid crystal in the mean field 
approximation. Being a low-dimensional model with nn 
interaction, the Id LL model does not exhibit any fi- 
nite temperature phase transition. This model has been 
solved exactly by Vuillermot and Romerio [38^ in 1973, 
using a group theoretical method. The results obtained 
in [111 are quoted below. The partition function Zn{K) 
for the TV-particle system is given by 



Zn{K) 



exp 



-NK 



(2) 



where K — 3/2T is a dimensionless quantity. D is the 
Dawson function [39j given by 



D{x) = exp(-.x^) / e" du 
Jo 



(3) 



The dimensionless internal energy Un{K), entropy 
Sn{K) and the specific heat Cn{K) are given by 

N 2 2 \ J ^ ' 



N 2 2 



N 



1 - ^3/2 



K- 



In K-^'^D{K^' 



D- 



(5) 
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2 



(6) 



We decided to choose this model to test the performance 
of WL algorithm using the suggested spin update scheme 
so that a comparison can be made with the exact results 
available for the model. 
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III. COMPUTATIONAL TECHNIQUES 

In the first part of tliis Section, we will describe the 
computational techniques used to determine the fluctu- 
ations in the energy histogram. In the later part of this 
Section, we will discuss the method for the new spin up- 
date scheme. 

Let us first explain the notations and symbols relevant 
to the present work. The saturation value of the energy 
histogram fluctuation in the /c"* iteration is represented 
by Let fk be the modification factor for the /c*'* 
iteration. One usually starts with a modification factor 
/ = /i > 1 and uses a sequence of decreasing fk's {k = 
1,2,3,---) defined in some manner. One Monte Carlo 
(MC) sweep is taken to be completed when the number of 
attempted single spin moves equals the number of spins in 
the system. The error in the DOS after the n*'' iteration 
is directly related to /3i for i > n, the saturation values 
of the fluctuations. In the WL algorithm the logarithm 
of the DOS after n iterations is given by 

n 

gr,{Ei) = J2Hk{E,)ln{fk) (7) 

where Hk{Ei) is the accumulated histogram count for 
the energy bin during the k^^ iteration. In order to 
get an idea of the fluctuations in the histogram and its 
growth with the number of MC sweeps, we subtract the 
minimum of the histogram count hi. which occurs in the 
histogram after the j*'' MC sweep has been completed 
during the /c*'* iteration, i.e., we consider the quantity 

HliEO ^ HiiEi) - hi (8) 

It may be noted that hj, does not refer to any particular 
bin and may occur in any of the visited bins. The quan- 
tity Hl{Ei) is now summed over all bins to give ^H^. 

i 

AiJ^ is thus a measure of the fluctuations which occurs 
in the j*'* MC sweep during A:*'' iteration and is a sort of 
average over all macrostates or bins. AiJ;^ fluctuates with 
j because of statistical errors and its mean value taken 
over j is nothing but /3fe. The error of the logarithm of 
the DOS, summed over all energy levels or bins, after the 
completion of n iterations is therefore given by [29j 

oo 

Vn= PkHfk) (10) 

k=n+l 

Eq. pH]) means that the error depends only on the fluc- 
tuations in histogram and the sequence of modification 
factors. When the values of fk are predetermined, the 
fluctuations in histogram, i.e., A77^, becomes the only 
determining factor for the error. For this reason the ob- 
servable Ai/^, defined by Eq. ([5]), is considered to be 



a good measure of the fiuctutations in histogram. How- 
ever, we point out that because of the summation over 
the index i in Eq. ([5]), the nature of the distribution 
of the errors over the energy bins is not reflected in the 
summed quantity A_ff^. What we get instead is an error 
which has been summed over all the energy bins. Since 
the predicted value of the error rjn is of the order of -\/ln /„ 
[27j . one expects that the histogram saturation value /3„, 
for the n*'* iteration, should be proportional to 1/ -^In /„ . 

A. Proposal for a novel spin update method 

Now we discuss the method to generate a subsequent 
less-correlated spin configuration. In the conventional 
spin update method for a continuous lattice spin model, 
the orientation of each spin s is stored in terms of the 
direction cosines (liThih)- To generate a new config- 
uration (microstate), a spin is selected at random and 
each direction cosine of it is updated as U ^ U + p * Xi 
for (j = 1,2,3) where the parameter "p" denotes the 
amplitude of the random angular displacements, chosen 
such that approximately half of the configurations are 
accepted and half rejected 0) and Xi IS cL random num- 
ber between —1 to 1. We have seen for a number of 
continuous lattice spin models that the results for the 
thermodynamic quantities become very sensitive to the 
value of the parameter "p" . "p" is generally taken such 
that p < 1 and the choice of "p" also depends on the 
systems we are working on. The reason for taking p < 1 
is that small values of "p" correspond to small changes 
in the direction of the spin, i.e., the energy cost of an 
attempted move will be small. However, this is not the 
only form of update, nor is it known whether this is the 
most efficient form. The thing is, there is quite a lot of 
fiexibility about the choice of the new state for the spins. 
A good discussion of it may be found in Ref . [ilj . 

In the present work, we propose a novel protocol to 
generate a less-correlated spin configuration in the fol- 
lowing manner. We take a random unit vector f and a 
spin update s — >■ s ' is defined as s ' = s — 2 (s • r) r where 
(s • r) is the dot product of s and f. This represents a 
refiection with respect to the hyperplane orthogonal to r 
and this is an idempotent operation. The idea came from 
Wolff . One may think of a linear transformation R{f) 
such that s ' = R{r)s. This linear transformation has the 
property 

Riff = 1 (11) 

i.e., idempotent and 

[i?(f)si] • [Rif}S2] = Si • S2 (12) 

i.e., the Hamiltonian ([Ij is invariant under global R 
transformations. This spin update method reduces the 
autocorrelation time to a considerable amount and con- 
sequently systematic error decreases. Moreover, defining 
a spin update in that way, the algorithm becomes free 
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from tuning any adjustable parameter even while simu- 
lating a lattice spin model with continuous energy spec- 
trum. This spin update method has resulted in efficient 
simulation of continuous lattice spin models with XY 
symmetry [isllii]. 

The energy of the Id LL model is a continuous variable 
and it can have any value between —L to L/2 where L 
is the system size. To discretize the system, we have 
chosen an energy range [—L, 0) and divided this energy 
range into a number of bins (macrostates) each having a 
width, say w. In the present work, the bin width is taken 
to be 0.2. 



IV. RESULTS AND DISCUSSIONS 

We have determined for the lattice model we have de- 
fined, the dependence of the quantity Ai/^ , given by Eq. 
([9]), on j, the number of MC sweeps for a given iteration 
denoted by k. For the purpose of testing the fluctuations 
in histogram, we have taken linear spin chains of length 
i = 80 and 160. Nearest neighbor interactions along 
with periodic boundary conditions were always used. The 
starting value of the modification factor In fi was taken 
to be 0.1 and the sequence ln/„+i = (ln/„)/10^/'' was 
chosen and for the purpose of determination of fiuctu- 
ations, the minimum In/ used was 10~^. Clearly, the 
chosen sequence of / is to ensure that it gets reduced by 
a factor of 10 after four iterations. We have determined 
the quantity AiJ^ defined by Eq. ^ at intervals of lO'^ 
MC sweeps and the maximum number of sweeps chosen 
for a given value of / is such that the saturation of the 
histogram is clearly evident. The system energy is always 
considered up to i? = 0. The lower limit of the energy 
for i = 80 is taken to be —78 and for L — 160, it is 
— 158, while the corresponding ground state energies are 
—80 and —160. Thus the visited energy range goes to a 
sufficiently low value to cover the entire range of interest, 
though the small cut near the ground state is necessary, 
as configurations near the minimum energy take a very 
long time to be sampled during the random walk. 

In Fig. [U we have plotted the fluctuations in the his- 
togram AJf]^ against the number of MC sweeps j for four 
values of the modification factor /. The plots shown are 
for L — 160 lattice and for In/ equal to 10~^, 10~^, 
10~* and 10~^. We did not go to values of In / less than 
10~^ as it takes a very large CPU time. Averages were 
taken over hundred independent simulations to improve 
the statistics and accuracy. Similar plots are also taken 
for the L = 80 lattice. It is evident from Fig. [T] that 
A_ff^ increases initially and then saturates and as / gets 
smaller, the saturation value as well as the number of 
MC sweeps necessary to reach the saturation (MCSsat) in- 
creases. Fig. [5] explicitly reveals this fact. The standard 
error calculated from the hundred independent simula- 
tions are also shown in Fig. [TJ In Fig. |31 we have plot- 
ted the logarithm of saturation value i.e., ln(/?fe) vs 
ln(ln /) for system sizes L = 80 and L = 160. From this 



figure, it is clear that 

/3fc(x(ln/)" (13) 

where the index a = -0.50133 ± 0.007 for i = 80 and 
a = -0.50844 ± 0.005 for L = 160 respectively. This 
is in agreement with the prediction of Zhou and Bhatt 
(27| . Certainly, this result is not new. It confirms the 
previous results that the values of the slope is generic to 
the WL algorithm, in this case, it is a continuous lattice 
spin model with spin dimensionality three. 

Now we present the results of various thermodynamic 
quantities obtained from the simulation. In Fig. Ul we 
have plotted the average energy per spin against tem- 
perature (T) for L = 220. The results have been com- 
pared with the exact values of this observable obtained 
from Ref. [ss']. The specific heat, calculated as fiuctu- 
ations of the energy, has been plotted against T in Fig. 
[5] for L = 220 and compared with the exact results. In 
the inset of Fig. [SJ the percentage error (e) in the Cy 
near the peak in comparison with the exact results is 
shown. Percentage error is a measure of how inaccurate 
(or accurate) a measurement is and is defined by the 

formula -^asured value-actual value ^ ^qqO/ g^act rCSultS 
actual value 

show that the specific heat peak is maximum at a tem- 
perature Tm^j. = 0.24 and from our simulation we obtain 
the temperature at which the peak of the specific heat is 
maximum is T^™ = 0.2351 for L = 220. This implies 
that the percentage error in temperature at which the 
peak of the specific heat is maximum is 2.04%. Fig. [6] 
shows the variation of entropy per particle for L = 220 
and the exact results are also shown in the same plot. 

The attention is now focused on the autocorrelation 
time of the simulation. The autocorrelation function for 
an observable 0{t) is given by 

Xit) = j dt\0{t') - {0)){0{t' + t)- (O)) (14) 

where 0{t) is the instantaneous value of the observable 
at time t and (O) is the average value. The integrand 
in the above equation actually measures the correlation 
between the fluctuation of the observables at two differ- 
ent times, one an interval t later than the other. So x(t) 
will take a nonzero value if on the average the fluctua- 
tions are correlated, otherwise it is zero. Thus when t is 
just a single MC step apart, we will have a large positive 
autocorrelation. For large x(i) will be zero and the 
measurements are totally uncorrelated. The autocorre- 
lation is expected to fall off exponentially at long times 
thus: 

X(<) ^ e-*/- (15) 

where r is a measure of autocorrelation time of our simu- 
lation. At time t = r, the autocorrelation function, which 
is a measure of the similarity of the two states, is only 
a factor of 1/e down from its maximum value at i = 0. 
We have estimated the autocorrelation time both for the 
simulations with the conventional spin update method 
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FIG. 1. (Color online) The histogram fluctuations Aff^ for the k*^ iteration are plotted against the MC sweeps j for the 
L = 160 lattice. The values of In/fe are indicated in the figures. The histograms are averaged over 100 independent simulations. 
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FIG. 2. (Color online) Logarithm of the MC sweeps required 
to reach the saturation is plotted against In/ for L = 160 
system size. The errorbars are shown in the figure. 



FIG. 3. (Color online) Logarithm of the saturation values of 
the histogram fiuctuations ln(/3fc) is plotted against ln(ln/) 
for L = 80 and L — 160 systems with the error bars. The 
slopes of the two linear fits are given in the text. 



and the proposed spin update method. The autocorrela- 
tion time is calculated following the method proposed by 
Madras and Sokal [1^. In the conventional spin update 
scheme, when we flip a single spin in each update, the 
total energy can only change by a small amount every 
time. In the proposed spin update scheme, the change in 
total energy is greater compared to that in the conven- 



tional scheme. As the WL algorithm does not require Ei 
to be close to Ef, but compares only il.{Ei) and il{Ef), 
the convergence becomes faster with the proposed scheme 
than with the conventional scheme. 

We have found that the autocorrelation time (r) ex- 
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FIG. 4. (Color online) The variation of the average energy 
per particle is plotted against T (solid line) for L = 220. 
Exact results, indicated by the filled circle, are also plotted in 
the same graph. The error bars are of the dimensions smaller 
than the symbols used for plotting. 



FIG. 6. (Color online) The variation of the entropy per par- 
ticle is plotted against T (solid line) for L = 220. The filled 
circle indicates the exact results. The error bars are of the 
dimensions smaller than the symbols used for plotting. 
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FIG. 5. (Color online) The specific heat is plotted against T 
(solid line) for L = 220. Exact results are indicated by the 
filled circle. The error bars are of the dimensions smaller than 
the symbols used for plotting. The percentage error in the Cv 
in comparison with the exact results is shown in the inset. 



hibits a power law scaling with system size, i.e., 

T cx (16) 

The scaling exponent (z) is determined from a linear fit 
of the plot Inr versus \a.L. The logarithm of the auto- 
correlation time for both the conventional and the pro- 
posed spin update scheme has been plotted against In L 
for In / = 0.01 in Fig. [T) The scaling exponent for the 
proposed spin update scheme (znew) is found to be Znew — 
1.36351 ± 0.024 while that for the conventional spin up- 
date scheme (zoid) is found to be Zoid = 1.57591 ± 0.013. 
The proposed spin update scheme significantly decreases 
the scaling exponent. We would like to point out that the 
autocorrelation time (r) increases rapidly as the modifi- 
cation factor (/) becomes smaller and for a larger system 
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FIG. 7. (Color online) Scaling of the autocorrelation time as 
a function of system size L for In / = 0.01 with the error bars 
shown. The scaling exponents are mentioned in the text. 



size, the calculation of r, specially for smaller /, becomes 
very much costly in terms of CPU time. The autocorre- 
lation time for a number of modification factors / for 
L ~ 200 for both the proposed and the conventional spin 
update schemes is listed in Table |T] and plotted in Fig. [51 



TABLE I. Autocorrelation time (in units of MC sweep) for 
different In / for L = 200. 



In/ 


Tnew 


Told 


1.0 


28954 


78966 


0.1 


31226 


96902 


0.01 


67642 


234502 


0.001 


173250 


638082 


0.0001 


246118 


1359694 
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FIG. 8. (Color online) Logarithm of the autocorrelation time 
plotted against In / for system size L — 200 for both the 
proposed and the conventional spin update method. 



V. CONCLUSIONS 

To summarize, we have tested the performance of the 
WL algorithm in a continuous lattice spin model, namely, 
the Id LL model which describes a nematic liquid crystal 
in the mean field approximation. The results obtained 
from our simulation are compared with the exact results 
available for this model. It has been observed that the 
results obtained tally accurately with the exact results. 
We focus on the fluctuations of histogram and replace 
the "flatness" criterion with that of minimum histogram. 
We have found that in this continuous lattice model, the 
fluctuations in the energy histogram, after an initial accu- 
mulation stage, saturates to a value that is proportional 
to 1/ \/\nf where / is the modification factor in the WL 
algorithm and confirm that this behavior is generic to 
the WL algorithm. We also present a novel method for 
spin update scheme to obtain a subsequent configuration 
which is less-correlated than the previous method. The 
proposed spin update scheme makes the WL "driver" 



to move from one sampling point to the next faster. 
As a result, the autocorrelation time between successive 
moves decreases and the convergence becomes faster. It 
may be noted that the WL algorithm only asks for the 
next sampling point (say X) with probability distribu- 
tion P{X) cx n{X)/n{X) where n{X) and n{X) are 
the exact and the estimated DOS respectively. A previ- 
ous study suggested that TV-fold way updates yields 
better performance in fiat-histogram sampling. However, 
Dayal et.al [1^ argued that the performance is limited by 
the added expense of the CPU time needed to implement 
the A^-fold way updates. The proposed method is simple 
to implement and has also the merit that it makes us free 
from tuning any adjustable parameter while simulating a 
continuous lattice spin model. Although the method has 
been applied to a liquid crystalline system in the present 
work, the method can, in general, be applied to any lat- 
tice spin model with continuous energy spectrum. This 
method has resulted in efficient simulation of continuous 
models with XY symmetry [i^ Finally, we stress 
that the focus in this paper is to test the performance of 
the WL algorithm in continuous lattice spin models with 
the proposed spin update scheme. We hope that this spin 
update method will be of general interest in the area of 
research in Monte Carlo simulations of continuous lattice 
spin models. 
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